home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / TPL60N14.ARJ / UNIT2.PAS < prev    next >
Pascal/Delphi Source File  |  1992-02-27  |  25KB  |  822 lines
  1. {$a+,n-,x-,s-,i-,r-,b-,v-}
  2.  
  3. unit Unit2;
  4. interface
  5.    uses mainvars;
  6.    procedure mile70170;
  7. implementation
  8.    procedure mile70170;
  9.    begin
  10.  
  11. {=============================================}
  12.    Milestone := 70;
  13. {=============================================}
  14.    writeln;
  15.    writeln ('Running test of square root(x).');
  16.    TestCondition (Failure, (Zero = sqrt (Zero))
  17.          and (- Zero = sqrt (- Zero))
  18.          and (One = sqrt (One)), ' Square root of 0.0, -0.0 or 1.0 wrong  '
  19.          );
  20.    MinSqrtError := Zero;
  21.    MaxSqrtError := Zero;
  22.    J := 0;
  23.    X := Radix;
  24.    OneUlp := U2;
  25.    SqrtXMinX (SeriousDefect);
  26.    X := BInverse;
  27.    OneUlp := BInverse * U1;
  28.    SqrtXMinX (SeriousDefect);
  29.    X := U1;
  30.    OneUlp := U1 * U1;
  31.    SqrtXMinX (SeriousDefect);
  32.    if J <> 0 then
  33.       begin
  34.       NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
  35.       Pause;
  36.       end;
  37.    writeln ('Testing if sqrt(X * X) = X for ', NoTrials, ' integers X.');
  38.    J := 0;
  39.    X := Two;
  40.    Y := Radix;
  41.    if (Radix <> One) then
  42.       repeat
  43.          X := Y;
  44.          Y := Radix * Y;
  45.       until (Y - X >= NoTrials);
  46.    OneUlp := X * U2;
  47.    I := 1;
  48.    Continue := true;
  49.    while (I <= NoTrials) and Continue do (* dgh: 10 --> NoTrials *)
  50.       begin
  51.       X := X + One;
  52.       SqrtXMinX (Defect);
  53.       if J > 0 then
  54.          begin
  55.          Continue := false;
  56.          NoErrors [Defect] := NoErrors [Defect] + 1;
  57.          end;
  58.       I := I + 1;
  59.       end;
  60.    writeln ('Test for Sqrt Monotonicity.');
  61.    I := - 1;
  62.    X := BMinusU2;
  63.    Y := Radix;
  64.    Z := Radix + Radix * U2;
  65.    NotMonot := false;
  66.    Monot := false;
  67.    while not (NotMonot or Monot) do
  68.       begin
  69.       I := I + 1;
  70.       X := sqrt (X);
  71.       Q := sqrt (Y);
  72.       Z := sqrt (Z);
  73.       if (X > Q) or (Q > Z) then
  74.          NotMonot := true
  75.       else
  76.          begin
  77.          Q := Int (Q + Half);
  78.          if (I > 0) or (Radix = Q * Q) then
  79.             Monot := true
  80.          else if I > 0 then
  81.             begin
  82.             if I > 1 then
  83.                Monot := true
  84.             else
  85.                begin
  86.                Y := Y * BInverse;
  87.                X := Y - U1;
  88.                Z := Y + U1;
  89.                end
  90.             end
  91.          else
  92.             begin
  93.             Y := Q;
  94.             X := Y - U2;
  95.             Z := Y + U2;
  96.             end
  97.          end
  98.       end;
  99.    if Monot then
  100.       writeln ('Sqrt has passed a test for Monotonicity.')
  101.    else
  102.       begin
  103.       NoErrors [Defect] := NoErrors [Defect] + 1;
  104.       writeln('DEFECT:  Sqrt(X) is non-monotonic for X near ', Y);
  105.       end;
  106. {=============================================}
  107.    Milestone := 80;
  108. {=============================================}
  109.    MinSqrtError := MinSqrtError + Half;
  110.    MaxSqrtError := MaxSqrtError - Half;
  111.    Y := (sqrt (One + U2) - One) / U2;
  112.    SqrtError := (Y - One) + U2 / Eight;
  113.    if SqrtError > MaxSqrtError then
  114.       MaxSqrtError := SqrtError;
  115.    SqrtError := Y + U2 / Eight;
  116.    if SqrtError < MinSqrtError then
  117.       MinSqrtError := SqrtError;
  118.    Y := ((sqrt (F9) - U2) - (One - U2)) / U1;
  119.    SqrtError := Y + U1 / Eight;
  120.    if SqrtError > MaxSqrtError then
  121.       MaxSqrtError := SqrtError;
  122.    SqrtError := (Y + One) + U1 / Eight;
  123.    if SqrtError < MinSqrtError then
  124.       MinSqrtError := SqrtError;
  125.    OneUlp := U2;
  126.    X := OneUlp;
  127.    for Index := 1 to 3 do
  128.       begin
  129.       Y := sqrt ((X + U1 + X) + F9);
  130.       Y := ((Y - U2) - ((One - U2) + X)) / OneUlp;
  131.       Z := ((U1 - X) + F9) * Half * X * X / OneUlp;
  132.       SqrtError := (Y + Half) + Z;
  133.       if SqrtError < MinSqrtError then
  134.          MinSqrtError := SqrtError;
  135.       SqrtError := (Y - Half) + Z;
  136.       if SqrtError > MaxSqrtError then
  137.          MaxSqrtError := SqrtError;
  138.       if ((Index = 1) or (Index = 3)) then
  139.          X := OneUlp * Sign (X) * Int (Eight / (Nine * sqrt (OneUlp)))
  140.       else
  141.          begin
  142.          OneUlp := U1;
  143.          X := - OneUlp;
  144.          end;
  145.       end;
  146. {=============================================}
  147.    Milestone := 85;
  148. {=============================================}
  149.    SquareRootWrong := false;
  150.    AnomolousArithmetic := false;
  151.    RSqrt := Other; (* ~dgh *)
  152.    if Radix <> One then
  153.       begin
  154.       writeln ('Testing whether sqrt is rounded or chopped: ');
  155.       D := Int (Half + Power (Radix, One + Precision - Int (Precision)))
  156.          ;
  157.    { ... = Radix^(1 + fract) if Precision = integer + fract. }
  158.       X := D / Radix;
  159.       Y := D / A1;
  160.       if (X <> Int (X)) or (Y <> Int (Y)) then
  161.          begin
  162.          AnomolousArithmetic := true;
  163.          end
  164.       else
  165.          begin
  166.          X := Zero;
  167.          Z2 := X;
  168.          Y := One;
  169.          Y2 := Y;
  170.          Z1 := Radix - One;
  171.          FourD := Four * D;
  172.          repeat
  173.             if Y2 > Z2 then
  174.                begin
  175.                Q := Radix;
  176.                Y1 := Y;
  177.                repeat
  178.                   X1 := abs (Q + Int (Half - Q / Y1) * Y1);
  179.                   Q := Y1;
  180.                   Y1 := X1;
  181.                until X1 <= Zero;
  182.                if Q <= One then
  183.                   begin
  184.                   Z2 := Y2;
  185.                   Z := Y;
  186.                   end;
  187.                end;
  188.             Y := Y + Two;
  189.             X := X + Eight;
  190.             Y2 := Y2 + X;
  191.             if Y2 >= FourD then
  192.                Y2 := Y2 - FourD;
  193.          until Y >= D;
  194.          X8 := FourD - Z2;
  195.          Q := (X8 + Z * Z) / FourD;
  196.          X8 := X8 / Eight;
  197.          if Q <> Int (Q) then
  198.             AnomolousArithmetic := true
  199.          else
  200.             begin
  201.             Break := false;
  202.             repeat
  203.                X := Z1 * Z;
  204.                X := X - Int (X / Radix) * Radix;
  205.                if X = One then
  206.                   Break := true
  207.                else
  208.                   Z1 := Z1 - One;
  209.             until Break or (Z1 <= 0);
  210.             if (Z1 <= 0) and (not Break) then
  211.                AnomolousArithmetic := true
  212.             else
  213.                begin
  214.                if Z1 > RadixD2 then
  215.                   Z1 := Z1 - Radix;
  216.                repeat
  217.                   NewD;
  218.                until U2 * D >= F9;
  219.                if D * Radix - D <> W - D then
  220.                   AnomolousArithmetic := true
  221.                else
  222.                   begin
  223.                   Z2 := D;
  224.                   I := 0;
  225.                   Y := D + (One + Z) * Half;
  226.                   X := D + Z + Q;
  227.                   SubRout3750;
  228.                   Y := D + (One - Z) * Half + D;
  229.                   X := D - Z + D;
  230.                   X := X + Q + X;
  231.                   SubRout3750;
  232.                   NewD;
  233.                   if D - Z2 <> W - Z2 then
  234.                      AnomolousArithmetic := true
  235.                   else
  236.                      begin
  237.                      Y := (D - Z2) + (Z2 + (One - Z) * Half);
  238.                      X := (D - Z2) + (Z2 - Z + Q);
  239.                      SubRout3750;
  240.                      Y := (One + Z) * Half;
  241.                      X := Q;
  242.                      SubRout3750;
  243.                      if I = 0 then
  244.                         AnomolousArithmetic := true;
  245.                      end
  246.                   end
  247.                end
  248.             end
  249.          end;
  250.       if (I = 0) or AnomolousArithmetic then
  251.          begin
  252.          NoErrors [Failure] := NoErrors [Failure] + 1;
  253.          writeln ('FAILURE:  Anomolous arithmetic with ',
  254.             'integer < Radix^Precision = ');
  255.          writeln (W, '  fails test whether sqrt rounds or chops.');
  256.          SquareRootWrong := true;
  257.          end
  258.       end;
  259.    if not AnomolousArithmetic then
  260.       begin
  261.       if not ((MinSqrtError < 0) or (MaxSqrtError > 0)) then
  262.          begin
  263.          RSqrt := Rounded;
  264.          writeln ('Square root appears to be correctly rounded.');
  265.          end
  266.       else if (MaxSqrtError + U2 > U2 - Half) or (MinSqrtError > Half)
  267.             or (MinSqrtError + Radix < Half) then
  268.          SquareRootWrong := true
  269.       else
  270.          begin
  271.          RSqrt := Chopped;
  272.          writeln ('Square root appears to be chopped.');
  273.          end;
  274.       end;
  275.    if SquareRootWrong then
  276.       begin
  277.       writeln ('Square root is neither chopped nor correctly rounded.');
  278.       writeln ('Observed errors run from ', MinSqrtError - Half);
  279.       writeln ('to ', Half + MaxSqrtError, ' ulps.');
  280.       TestCondition (SeriousDefect, MaxSqrtError - MinSqrtError < Radix
  281.             * Radix, 'sqrt gets too many last digits wrong.   ');
  282.       end;
  283. {=============================================}
  284.    Milestone := 90;
  285. {=============================================}
  286.    Pause;
  287.    (* fix from dgh: Wichman had effectively commented out here to Milestone 110 *)
  288.       writeln ('Testing powers Z^i for small integers Z and i.');
  289.       N := 0;
  290.    { ... test power of zero. }
  291.       I := 0;
  292.       Z := - Zero;
  293.       M := 3;
  294.       Break := false;
  295.       repeat
  296.          X := One;
  297.          SubRout3980;
  298.          if I <= 10 then
  299.             begin
  300.             I := 1023;
  301.             SubRout3980;
  302.             end;
  303.          if Z = MinusOne then
  304.             Break := true
  305.          else
  306.             begin
  307.             Z := MinusOne;
  308.          { .. if(-1)^N is invalid, replace MinusOne by One. }
  309.             I := - 4;
  310.             end;
  311.       until Break;
  312.       PrintIfNPositive;
  313.       N1 := N;
  314.       N := 0;
  315.       Z := A1;
  316.       M := Int (Two * ln (W) / ln (A1));
  317.       Break := false;
  318.       repeat
  319.          X := Z;
  320.          I := 1;
  321.          SubRout3980;
  322.          if Z = AInverse then
  323.             Break := true
  324.          else
  325.             Z := AInverse;
  326.       until Break;
  327.    {=============================================}
  328.       Milestone := 100;
  329.    {=============================================}
  330.    {  Power of Radix have been tested, }
  331.    {         next try a few primes     }
  332.       M := NoTrials;
  333.       Z := Three;
  334.       repeat
  335.          X := Z;
  336.          I := 1;
  337.          SubRout3980;
  338.          repeat
  339.             Z := Z + Two;
  340.          until (Three * Int (Z / Three) <> Z);
  341.       until (Z >= Eight * Three);
  342.       if N > 0 then
  343.          begin
  344.          writeln ('Error like this may invalidate financial ');
  345.          writeln ('calculations involving interest rates.');
  346.          end;
  347.       PrintIfNPositive;
  348.       N := N + N1;
  349.       if N = 0 then
  350.          writeln ('... no discrepancies found.');
  351.       writeln;
  352.       if N > 0 then
  353.          Pause;
  354. {=============================================}
  355.    Milestone := 110;
  356. {=============================================}
  357.    writeln ('Seeking Underflow thresholds UnderflowThreshold and E0');
  358.    D := U1;
  359.    if (Precision <> Int (Precision)) then
  360.       begin
  361.       D := BInverse;
  362.       X := Precision;
  363.       repeat
  364.          D := D * BInverse;
  365.          X := X - One;
  366.       until X <= Zero;
  367.       end;
  368.    Y := One;
  369.    Z := D;
  370. { ... D is power of 1/Radix < 1. }
  371.    repeat
  372.       C := Y;
  373.       Y := Z;
  374.       Z := Y * Y;
  375.    until not ((Y > Z) and (Z + Z > Z));
  376.    Y := C;
  377.    Z := Y * D;
  378.    repeat
  379.       C := Y;
  380.       Y := Z;
  381.       Z := Y * D;
  382.    until not ((Y > Z) and (Z + Z > Z));
  383.    if Radix < Two then
  384.       HInverse := Two
  385.    else
  386.       HInverse := Radix;
  387.    H := One / HInverse;
  388. { ... 1/HInverse = H = Min(1/Radix, 1/2) }
  389.    CInverse := One / C;
  390.    E0 := C;
  391.    Z := E0 * H;
  392. { ...1/Radix^(BIG integer) << 1 << CInverse = 1/C }
  393.    repeat
  394.       Y := E0;
  395.       E0 := Z;
  396.       Z := E0 * H;
  397.    until not ((E0 > Z) and (Z + Z > Z));
  398.    UnderflowThreshold := E0;
  399.    E1 := Zero;
  400.    Q := Zero;
  401.    E9 := U2;
  402.    S := One + E9;
  403.    D := C * S;
  404.  
  405.    if D <= C then
  406.       begin
  407.       E9 := Radix * U2;
  408.       S := One + E9;
  409.       D := C * S;
  410.       if D <= C then
  411.          begin
  412.          writeln ('FAILURE:  multiplication gets too many last digits wrong.');
  413.          NoErrors [Failure] := NoErrors [Failure] + 1;
  414.          Underflow := E0;
  415.          Y1 := Zero;
  416.          PseudoZero := Z;
  417.          Pause;
  418.          end
  419.       end
  420.    else
  421.       begin
  422.       Underflow := D;
  423.       PseudoZero := Underflow * H;
  424.       UnderflowThreshold := Zero;
  425.       repeat
  426.          Y1 := Underflow;
  427.          Underflow := PseudoZero;
  428.          if E1 + E1 <= E1 then
  429.             begin
  430.             Y2 := Underflow * HInverse;
  431.             E1 := abs (Y1 - Y2);
  432.             Q := Y1;
  433.             if (UnderflowThreshold = Zero) and (Y1 <> Y2) then
  434.                UnderflowThreshold := Y1;
  435.             end;
  436.          PseudoZero := PseudoZero * H;
  437.       until not ((Underflow > PseudoZero)
  438.             and (PseudoZero + PseudoZero > PseudoZero));
  439.       end;
  440. { Comment line 4530 .. 4560 }
  441.    if PseudoZero <> Zero then
  442.       begin
  443.       writeln;
  444.       Z := PseudoZero;
  445.    { ... Test PseudoZero for "phoney-zero" violating }
  446.    { ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
  447.          ... }
  448.       if PseudoZero <= 0 then
  449.          begin
  450.          NoErrors [Failure] := NoErrors [Failure] + 1;
  451.          writeln ('FAILURE:  Positive expressions can underflow to an ');
  452.          writeln ('allegedly negative value PseudoZero that prints out as:');
  453.          writeln (PseudoZero);
  454.          X := - PseudoZero;
  455.          if X <= 0 then
  456.             begin
  457.             writeln ('But -PseudoZero, which should then be positive, isn''t;');
  458.             writeln (' it prints out as: ', X);
  459.             end
  460.          end
  461.       else
  462.          begin
  463.          NoErrors [Flaw] := NoErrors [Flaw] + 1;
  464.          writeln ('FLAW:  Underflow can stick at an allegedly positive');
  465.          writeln ('value PseudoZero that prints out as: ', PseudoZero);
  466.          end;
  467.       TestPartialUnderflow;
  468.       end;
  469. {=============================================}
  470.    Milestone := 120;
  471. {=============================================}
  472.    if (CInverse * Y > CInverse * Y1) then
  473.       begin
  474.       S := H * S;
  475.       E0 := Underflow;
  476.       end;
  477.    if not ((E1 = 0) or (E1 = E0)) then
  478.       begin
  479.       NoErrors [Defect] := NoErrors [Defect] + 1;
  480.       write ('DEFECT:  ');
  481.       if E1 < E0 then
  482.          begin
  483.          writeln ('Products underflow at a higher');
  484.          writeln ('threshold than differences.');
  485.          if PseudoZero = Zero then
  486.             E0 := E1;
  487.          end
  488.       else
  489.          begin
  490.          writeln ('Difference underflows at a higher');
  491.          writeln ('threshold than products.');
  492.          end;
  493.       end;
  494.    writeln ('Smallest strictly positive number found is E0 =', E0);
  495.    Z := E0;
  496.    TestPartialUnderflow;
  497.    Underflow := E0;
  498.    if N = 1 then
  499.       Underflow := Y;
  500.    I := 4;
  501.    if E1 = Zero then
  502.       I := 3;
  503.    if UnderflowThreshold = Zero then
  504.       I := I - 2;
  505.    UnderflowNotGradual := true;
  506.    case I of
  507.       1:
  508.          begin
  509.          UnderflowThreshold := Underflow;
  510.          if (CInverse * Q) <> ((CInverse * Y) * S) then
  511.             begin
  512.             NoErrors [Failure] := NoErrors [Failure] + 1;
  513.             UnderflowThreshold := Y;
  514.             writeln ('FAILURE:  Either accuracy deteriorates as numbers');
  515.             writeln ('approach a threshold UnderflowThreshold = ',
  516.                   UnderflowThreshold);
  517.             writeln ('coming down from ', C, ' or else multiplication');
  518.             writeln ('gets too many last digits wrong.');
  519.             end;
  520.          Pause;
  521.          end;
  522.       2:
  523.          begin
  524.          NoErrors [Failure] := NoErrors [Failure] + 1;
  525.          writeln ('FAILURE:  Underflow confuses Comparison, which alleges that');
  526.          writeln ('Q = Y while denying that |Q - Y| = 0 ; these values');
  527.          write ('print out as Q = ', Q, ' Y = ', Y2, ' |Q - Y| = ');
  528.          writeln (abs (Q - Y2));
  529.          UnderflowThreshold := Q;
  530.          end;
  531.       3:
  532.          begin
  533.          X := X;
  534.       {FOR PRETTYPRINTER, IS DUMMY}
  535.          end;
  536.       4:
  537.          if (Q = UnderflowThreshold) and (E1 = E0)
  538.                and (abs ( UnderflowThreshold - E1 / E9) <= E1) then
  539.             begin
  540.             UnderflowNotGradual := false;
  541.             writeln ('Underflow is gradual; it incurs Absolute Error =')
  542.                ;
  543.             writeln ('(roundoff in UnderflowThreshold) < E0.');
  544.             Y := E0 * CInverse;
  545.             Y := Y * (OneAndHalf + U2);
  546.             X := CInverse * (One + U2);
  547.             Y := Y / X;
  548.             IEEE := (Y = E0);
  549.             end;
  550.       end;
  551.    if UnderflowNotGradual then
  552.       begin
  553.       writeln;
  554.       R := sqrt (Underflow / UnderflowThreshold);
  555.       if R <= H then
  556.          begin
  557.          Z := R * UnderflowThreshold;
  558.          X := Z * (One + R * H * (One + H));
  559.          end
  560.       else
  561.          begin
  562.          Z := UnderflowThreshold;
  563.          X := Z * (One + H * H * (One + H));
  564.          end;
  565.       if not ((X = Z) or (X - Z <> 0)) then
  566.          begin
  567.          NoErrors [Flaw] := NoErrors [Flaw] + 1;
  568.          writeln ('FLAW:  X = ', X, ' is unequal to Z = ', Z);
  569.          write ('yet X - Z yields ');
  570.          Z9 := X - Z;
  571.          writeln (Z9, '.  Should this NOT signal Underflow,');
  572.          writeln ('this is a SERIOUS DEFECT that causes');
  573.          writeln ('confusion when innocent statements like');
  574.          write ('    if X = Z then  ...  else');
  575.          writeln ('  ...  (f(X) - f(Z)) / (X - Z) ...');
  576.          writeln ('encounter Division by Zero although actually');
  577.          writeln ('X / Z = 1 + ', (X / Z - Half) - Half);
  578.          end;
  579.       end;
  580.    writeln ('The Underflow threshold is ', UnderflowThreshold,
  581.          ' below which');
  582.    writeln ('calculation may suffer larger Relative error then',
  583.          ' merely roundoff.');
  584.    Y2 := U1 * U1;
  585.    Y := Y2 * Y2;
  586.    Y2 := Y * U1;
  587.    if Y2 <= UnderflowThreshold then
  588.       begin
  589.       if Y > E0 then
  590.          begin
  591.          NoErrors [Defect] := NoErrors [Defect] + 1;
  592.          I := 5;
  593.          end
  594.       else
  595.          begin
  596.          NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
  597.          write ('SERIOUS ');
  598.          I := 4;
  599.          end;
  600.       writeln ('DEFECT: Range is too narrow; U1^', I, ' Underflows.');
  601.       end;
  602. {=============================================}
  603.    Milestone := 130;
  604. {=============================================}
  605. { SKIP
  606.    Y := - Int (Half - TwoForty * ln (UnderflowThreshold) / ln (HInverse)
  607. ) / TwoForty;
  608.    Y2 := Y + Y;
  609.    writeln ('Since underflow occurs below the threshold');
  610.    writeln ('UnderflowThreshold = ( ', HInverse, ' ) ^ ( ', Y, ') only u
  611. nderflow');
  612.    writeln ('should afflict the expression ( ', HInverse, ' ) ^ ( ', Y2,
  613.          ' )');
  614.    write ('actually calculating yields: ');
  615.    V9 := Power (HInverse, Y2);
  616.    writeln (V9);
  617.    if not ((V9 >= 0) and (V9 <= (Radix + Radix + E9) * UnderflowThreshol
  618. d)) then
  619.       begin
  620.       NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
  621.       writeln ('SERIOUS DEFECT:  this is not between 0 and');
  622.       writeln (' UnderflowThreshold = ', UnderflowThreshold);
  623.       end
  624.    else if not (V9 > UnderflowThreshold * (One + E9)) then
  625.       writeln ('this computed value is O.K.')
  626.    else
  627.       begin
  628.       NoErrors [Defect] := NoErrors [Defect] + 1;
  629.       writeln ('DEFECT:  This is not between 0 and UnderflowThreshold = ',
  630.                UnderflowThreshold);
  631.       end;
  632. end SKIP }
  633. {=============================================}
  634.    Milestone := 140;
  635. {=============================================}
  636.    writeln;
  637. { ...calculate Exp2 = exp(2) = 7.389056099... }
  638.    X := 0;
  639.    I := 2;
  640.    Y := Two * Three;
  641.    Q := 0;
  642.    N := 0;
  643.    repeat
  644.       Z := X;
  645.       I := I + 1;
  646.       Y := Y / (I + I);
  647.       R := Y + Q;
  648.       X := Z + R;
  649.       Q := (Z - X) + R;
  650.    until X <= Z;
  651.    Z := (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
  652.    X := Z * Z;
  653.    Exp2 := X * X;
  654.    X := F9;
  655.    Y := X - U1;
  656.    write ('Testing X^((X + 1) / (X - 1)) vs. exp(2) = ');
  657.    writeln (Exp2, ' as X -> 1.');
  658.    Break := false;
  659.    I := 1;
  660.    while (not Break) do
  661.       begin
  662.       Z := X - BInverse;
  663.       Z := (X + One) / (Z - (One - BInverse));
  664.       Q := Power (X, Z) - Exp2;
  665.       if abs (Q) > TwoForty * U2 then
  666.          begin
  667.          Break := true;
  668.          N := 1;
  669.          NoErrors [Defect] := NoErrors [Defect] + 1;
  670.          writeln ('DEFECT:  Calculated ', Power(X,Z), ' for');
  671.          writeln (' (1 + (', (X - BInverse) - (One - BInverse),
  672.               ')) ^ (', Z, ')');
  673.          writeln (' which differs from the correct value by Q =', Q);
  674.          writeln (' This much error may spoil financial',
  675.               ' calculations involving');
  676.          writeln (' tiny interest rates.');
  677.          end
  678.       else
  679.          begin
  680.          Z := (Y - X) * Two + Y;
  681.          X := Y;
  682.          Y := Z;
  683.          Z := One + (X - F9) * (X - F9);
  684.          if ((Z > One) and (I < NoTrials)) then
  685.             I := I + 1
  686.          else if X > One then
  687.             begin
  688.             if N = 0 then
  689.                writeln ('Accuracy seems adequate.');
  690.             Break := true;
  691.             end
  692.          else
  693.             begin
  694.             X := One + U2;
  695.             Y := U2 + U2 + X;
  696.             I := 1;
  697.             end
  698.          end
  699.       end;
  700.  
  701. {=============================================}
  702.    Milestone := 150;
  703.    {=============================================}
  704.    writeln ('Testing powers Z^Q at four nearly extreme values.');
  705.    N := 0;
  706.    Z := A1;
  707.    Q := Int (Half - ln (C) / ln (A1));
  708.    Break := false;
  709.    repeat
  710.       X := CInverse;
  711.       Y := Power (Z, Q);
  712.       DoesYequalX;
  713.       Q := - Q;
  714.       X := C;
  715.       Y := Power (Z, Q);
  716.       DoesYequalX;
  717.       if Z < One then
  718.          Break := true
  719.       else
  720.          Z := AInverse;
  721.    until Break;
  722.    PrintIfNPositive;
  723.    if N = 0 then
  724.       writeln (' ... no discrepancies found.');
  725.    writeln;
  726.    if N > 0 then
  727.       Pause;
  728. {=============================================}
  729.    Milestone := 160;
  730. {=============================================}
  731.    Pause;
  732.    writeln ('Searching for Overflow threshold:  ',
  733.         'This may generate an error.');
  734.    writeln ('Try a few values for N, starting with a large one,');
  735.    writeln ('and take the one that just does not stop the machine.');
  736.    writeln ('Did you find the correct value for N yet?');
  737.    writeln ('Hint: the number for Turbo-Pascal is 45');
  738.    readln (input);
  739.    while not eoln (input) do
  740.       read (input, ch);
  741.    Break := false;
  742.    repeat
  743.       writeln ('N = ');
  744.       readln (input);
  745.       while not eoln (input) do
  746.          read (input, NoTimes);
  747.       Y := - CInverse;
  748.       V9 := HInverse * Y;
  749.       Index := 1;
  750.       I := 0;
  751.       repeat
  752.          V := Y;
  753.          Y := V9;
  754.          V9 := HInverse * Y;
  755.          if V9 >= Y then begin I := 1; ch := 'Y'; Index := NoTimes end;
  756.          Index := Index + 1;
  757.          until Index > NoTimes;
  758.       if I = 0 then Y := V9;
  759.  
  760.       if (ch = 'N') or (ch = 'n') then
  761.          writeln ('N seems not large enough, try again.')
  762.       else
  763.          begin
  764.          writeln ('O.K.');
  765.          Break := true;
  766.          end;
  767.    until Break;
  768.    Z := V9;
  769.    writeln ('Can "Z = -Y" overflow?  Trying it on Y = ', Y);
  770.    V9 := - Y;
  771.    V0 := V9;
  772.    if (V - Y = V + V0) then
  773.       writeln ('Seems O.K.')
  774.    else
  775.       begin
  776.       NoErrors [Flaw] := NoErrors [Flaw] + 1;
  777.       writeln ('Finds a FLAW:  -(-Y) differs from Y.');
  778.       end;
  779.    if Z <> Y then
  780.       begin
  781.       NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
  782.       writeln ('SERIOUS DEFECT:');
  783.       writeln ('  overflow past ', Y, ' shrinks to ', Z);
  784.       end;
  785.    if (I = 1) then
  786.       begin
  787.       Y := V * (HInverse * U2 - HInverse);
  788.       Z := Y + ((One - HInverse) * U2) * V;
  789.       if Z < V0 then
  790.          Y := Z;
  791.       if Y < V0 then
  792.          V := Y;
  793.       if V0 - V < V0 then
  794.          V := V0;
  795.       end
  796.    else begin
  797.       V := Y * (HInverse * U2 - HInverse);
  798.       V := V + ((One - HInverse) * U2) * Y;
  799.       end;
  800.    writeln ('Overflow threshold is V = ', V);
  801.    if (I = 1) then writeln ('Overflow saturates at V0 = ', V0)
  802.       else begin writeln('There is no saturation value because');
  803.           writeln('the system traps on overflow.') end;
  804.    write ('No Overflow should get signaled for V * 1 = ');
  805.    V9 := V * One;
  806.    writeln (V9);
  807.    write ('                            nor for V / 1 = ');
  808.    V9 := V / One;
  809.    writeln (V9);
  810.    writeln ('Any overflow signal separating this * from the one');
  811.    writeln ('above is a DEFECT.');
  812. {=============================================}
  813.    Milestone := 170;
  814. {=============================================}
  815.    TestCondition (Failure, (- V < V) and (- V0 < V0)
  816.          and (- UnderflowThreshold < V)
  817.          and (UnderflowThreshold < V),
  818.            'Comparisons are confused by Overflow.   ');
  819.  
  820.    end;
  821. end.
  822.